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Consider two simultaneous first order differential equations x'(t)=F(x,y,t), y'(t) = 
G(x,y,t). Runge-Kutta type integration methods are developed which allow differenl inte- 
gration steps to be used for these equations. These methods retain the desirable properties 
of Runge-Kutta methods, namely the self -starting property and ease of change of integration 
step. Two different approaches are considered and extensive experimental work is reported 
upon. Experiments are done both in situations where these methods are advantageous 
and where they are not. It is seen that these methods are more efficient than the normal 
Runge-Kutta methods if they are at all applicable and in ideal situations they give the same 
accuracy with 90 percent less computation. These methods are applicable to six degree 
of freedom missile simulations, for which they were developed. 

1. Introduction 
Consider two simultaneous first order differential equations: 

^=F(x,y,t) x(to)=x (1.1) 

^=G(x,y,t) y(t )=y (1.2) 

where y(t) does not depend strongly on x(t) or varies mueh more rapidly than x(t). In a 
normal numerical integration method for these equations, the integration step h must be chosen 
small enough to adequately integrate both (1.1) and (1.2). In this paper Runge-Kutta type 
methods are described which allow different integration steps to be used for these equations. 
These methods retain the desirable properties of Runge-Kutta methods, namely the self- 
starting property and the ability to change the integration step easily. 2 

The problem is defined in detail and two different approches to the development of the 
formulas are given in section 2. The analysis for third order integration formulas is given in 
sections 3, 4, 5, and 6. In section 7 the results are stated without derivation for fourth order 
integration. 

Three systems of differential equations have been solved using these formulas witli varying 
values of the parameters. The first equation is of the type suited for these formulas and they 
result in a considerable saving in computation. The results are discussed in detail in section 9. 
They point out that the second approach gives formulas which are considerably more accurate 
than the first approach — although one would not expect this beforehand. The second equation 
is of a type not suited for these formulas. The results are discussed in section 10. The third 
equation is oi the type suited for these formulas except that a very high frequency low-amplitude 
oscillation has been added to x(t). The experimental results are somewhat erratic. 

In the final section a detailed discussion is given for the situations where these formulas 
are most useful and also comparison is made of their efficiency with that of the usual methods. 
It is seen that these formulas are more efficient than the normal Runge-Kutta methods if they 
are at all applicable and that in ideal situations they may give the same accuracy as normal 
Runge-Kutta with 90 percent less computation. One area of application is to six degree of 
freedom missile simulations, for which these formulas were originally derived. 

i This work was done while the author was at Autonetics, Inc. and at the National Bureau of Standards as an NRC-NBS Research Associate. 
2 For an account of the basic properties of Runge-Kutta see F. B. Hildebrand, Introduction to Numerical Analysis (McGraw-Hill Book Co., 
Inc., New York, N.Y., 1956). 
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2. Problem Definition 



h denotes the integration step for (1.2), and Kh, where K is an integer, is the integration 
step for (1.1). Let t n denote t -\-nh and let x n and y n denote the numerical solutions of (1.1) 
and (1.2), respectively at t=t n . Round-off error is not considered in this paper. 

Consider two £-axes, 

<- Kh > 



kim+DK 



t ( 



m+2)K 



IsmK+l 



L(m+1)K-1 



l<(m+l)K+l 



U 



m+2)K 



the first for (1.1) and the second for (1.2). It is now desired to obtain Runge-Kutta type 
integration formulas that integrate (1.1) in steps of Kh and (1.2) in steps of h. 

(1.1) can be integrated from t mK to t (m +nK by a normal Runge-Kutta method. For a 
third order method the pertinent equations are: 

k =KhF(x n , y n , t n ) 

k 1 =KkF(x n +y 1 k 0j yn+yiK t n +yiKh) 

k 2 =KhF(x n +y s k 1 + (y 2 — y 3 )k , y n +y3h 1 + (y 2 —y 3 )h , t n +y 2 Kh) 

h =KhG(x n , y n , t n ) Y (2.1) 

h l =KhG(x n +y 1 k , Vn+yiK t n +y x Kh) 

h 2 =KhG(x n +y^k 1 + (y 2 —y i )k , y n +yzh l J r (72—73)^0, t n +y 2 Kh) 

X( m +i)K = %mK I A)^0T&^1 I P2^2- 

The integration parameters y u y 2 , 73, /3 , ft, and /3 2 may be those of any third order Runge- 
Kutta method. Note that h 2 need not be computed for this integration. 

The main difficulty in integrating (1.2) is to obtain values of x(t) at the integration points 
between t mK and < (m+1)x . A natural way to obtain these values is to extrapolate x(t) from t mK 
The Runge-Kutta method is itself an extrapolation process and one extrapolation has been 
made in the integration of (1.1) from t mK to t im+1)K . The values of k , k u and k 2 from (2.1) 
may also be used to extrapolate x(t) from t mK to the intermediate points. Let x mK+j denote the 
extrapolated value of x(t) at t mK+j , 1 <j<K-l. Then, for the appropriate coefficients \i(j), 

3 

Other estimates of x(t) are needed and they are obtained in the same manner. 

For a third order method the equations for integrating (1.2) from t mK +j to t mK + w are: 

^o (j) =hG (x mK+ j,y mK +j,t m K+j) 
d 1 (j)=hG(x mK+j +^\ i (j)k i -4jymK+j+^ido,t m K+3+^ih) 

T 

d 2 (i) =hG ( x mK+j +it] Xj (j) ki- 7 ,y mK +j+^di+ (^2— Ms) d ,t mK +j+H2fi>) 

y m K+j+i=ymK+j+cxAU)+ a idi (i)+«2^2 (j). 
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Two basic approaches are considered for obtaining the integration parameters p h /x 2 , /x 3 , 
a , ax and a 2 and the extrapolation parameters X x (j), . . . , X 9 (j). The first is to consider the 
truncation error of the resulting integration formula for y(t). The parameters can be chosen 
so that the method is third order. This is done in sections 3, 4, and 5. The second approach 
is to determine the extrapolation parameters so as to make the extrapolation as accurate as 
possible and then to determine the integration parameters independently. This is dour in 
section 6. 

3. Preliminary Computations 

For simplicity 2*tlK(j)ki- n will be denoted by 2 W . Likewise, 



A»(j) 

Tn(j) 

4.0") 

<t>n{j) 



7lX»+l&')+72Xn+20") 

y 2 iK+i(j)+y\K+2(j) 
7173X7H-2O') 



These expressions will occur quite often. The argument j will be omitted if it leads to no 
confusion. 

In order to consider the truncation error of the integration formula, various terms will be 
expanded in Taylor's series. There will be some derivatives evaluated at (x mK ,y mK ,t mK ) and 
some at (x mK +j,y m K+j,t mK +j)' Expansions will be obtained in this section that will allow a 
comparison of such terms. A Taylor's series expansion of S w will also be obtained. 

The notation: F p =bF/bp, G p =bG/bp, F PQ =b 2 F/bpdq, G pq =b 2 G/c>pi>q) P,q=x,y,t will 
be used. The convention that F denotes F(x mK ,y mK ,t mK ) and that G denotes G(x mK +j,y mK+j) 
tmK+j) will be adopted. The same convention will hold for derivatives of F and G. Further- 
more F X) G X) F 2 ,G 2 will denote dF/dt, dG/dt, 

F xx F 2 +2F xv FG+F yy G*+2F xt F+2F vt G+F tt , 

G xx F 2 +2G xy FG+G yv G*+2G xt F+2G yt G+G tt , 

respectively. 

Then it is seen that 

ft \X m K+j,ymK+j)tmK+j) 

= ft ~T ft x (XmK+j X mK ) -\~fty (y m K+j ymKj T^/jHll^M \X mK+ j X m Kj 
-\-2F xy (X mK +j X mK ) (y m K+j VrnKj ~rft yy\VmK+i VrnKj ~\~F tt (jh) 

+2F xt (x mK+j —x mK )(jh)+2F yt (y mK + j —y mK )(jh)} + . . . 

G(x m K,y m K,tmKj = G (xx{XmK+j x mK ) (j y (y rnK +j y m Kj G t {jh)-\- . . . 

(KhY r (3-1) 

S n = A n KhF+ (Kh) 2 I\[F x F+F y G(x Km y mK ,t mK ) +F t ] + K -f- (A w ) [F XX F 2 

-\-2r xy t (j(X m K,ymK)tmKJ ~T ft yyG \%mK,y mKitrnK) \2ft x tft ~l~2r yt G(X mK ,y mK ,t mK ) Tf « t] 

+ (Kh) 3 cf> n [F x (F x F+F y G(x mK ,y mK) t mK ) +F t )+F y (G x (x mK ,y mK ,t mK )F 

^(7 y (X m K,ymK,tmK)G(X mK ,y mK ,t mK )^G t (X m K . • . 

(»=l,4,7) 
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Let y n (t) denote the solution of (1.2) which assumes the value y n at t n . Then 

VmK+j VmK = ymK+j ymK+jV"mK) ~TymK+j\tmK) VmK' 

Since this will be a third order integration method, 

ymK+j(tmK)—ymK = ah*+bh 5 + .... 



Now 



and hence 



(jh) 2 
y m K + j-ymK+j(t n )=Gjh+^^[G x F(x mK+j ,y mK+j) t mK . hj )+G y G+G . . 



(3.2) 



07*/) 2 
y m K+j—ymK=Gjh-\ — k~ [G x F(x mK+j) y mK+ j } t mK +j)-\-GyG-\rGt]+ . . . 

With equation (3.2) and by repeated substitution, eqs (3.1) assume the form: 
F(x mK+j , y mK+j , t mK+j )=F+A 1 KhF x F+jh(F y G+F t ) 

+^[F xx F 2 K 2 A\+2F xy FGKA x j+F tt ^^ . . . (3.3) 



G(x mK ,y mK ,t mK ) = G-G x FKA 1 h-jh(G y G+G t ) + . . . 
? n =FKAji+(Kh) 2 Y n F x + { ^ A n F 2 +(Kh)*<l> n (F x F l +F 1 1 ) 



(3.4) 



(7i=l,4,7). (3.5) 



-(Kh) 2 T n [G x FKA 1 h+G y Gjh+G t jh} + . . . 

4. The First Approach 

The difference y mK+j (t mK+j+1 )—y mK+j+1 will be expanded in a Taylor's series and the coef- 
ficients of all terms of third order or less in h will be equated to zero. The resulting equations 
will be used to determine the integration and extrapolation parameters. Now 

ymK+j{tmK+j+l) ymK+j+l = ymK+j(tmK+j+l) ymK+j <*(A) #1^1 «2^2- 

With eqs (3.3), (3.4), and (3.5) it is seen that: 

y mK+j(tmK+j+i) = ymK+j + Gh+-h 2 [G x F(x mK+j ,y mK+j ,t mK + j ) -\-G y G-\- G t ] 

-j--/i [G tt -\-2(j ty G-\-GyyG -\-G xx r (x m K+j,y m K+jfmK+j) 
-\-2G xt r (x mK +j,y mK j r j,t mK j r j)-\r2G xy Gr (x mK+j ,y mK+j ,t mK+j )] 



d (j)=hG, 



+^h z [Gy(G x F(x mK+h y mK ^j,t mK+j ) + G y G+ G t ) 

-\- G x {F x (x mK +jjy mK+ j,t mK +j)F(x mK+ j,y mK+j ,t mK+ j) 

-f- GF y{x mK+ j ) y mK +j ) tmK+j) -\- F t (x mK +j,y mK+j ,t mK+ j) } ] + 

=y mK+} +Gh+\b?G l +\vG i 
+\b?{G v G x + Q^+^G^FFJKk, + F v Gj+F t j) + . 
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(4.1) 



(4.2) 



d 1 (j)=hG+h(G x 2 i +G y G^ 1 h+G tf x 1 h) 

+ I^[^21+2^(?S 4/ . 1 A+2^ ( 2 4/ z i A+^ ; , / ^( Mi / 1 )M-2^^(m i A) 2 +(?„(mi/i) 2 ] 
^hO+WiGsFKAi+GvGvu+GM) 

+ W[G xx F 2 K 2 Al+2G xy GFKA ifil +2G xl FKA 4 n 1 + G vy G* IJ .i 
+2G vt Gnl+G ltl 4]+VK>r i G x F l + ... 

d 2 (j) =hG+h(G x 2 7 +G y G,x 2 h+ Guijt) 

+ n 3 VG y (G x Z i +G v Gn 1 h+G tf ji 1 h)+ . . . 
=hG+h 2 [G x FKA 7 +G y G^ + G tt i 2 ] 
+ W[O xx F 2 K 2 Al+2G xy GFKA 7 n 2 +2G xl FKA^ 2 +G„G^l+2G yl G^ + G llf xl] 
+ n i h*G y (G x FKA i +G y Gm + G t n l )+h*K*r 7 G x F 1 + . . . 



(4.3) 



(4.4) 



The following system of equations results when the coefficients of the various terms are 
set equal to zero: 

aiMi+a2M2=o 



(4.5) 


ao+ai4-a 2 =l 


I 


aiju 1 4-a 2 /i2=2 


(4.6) 


a 1 A t K+a 2 A 7 K=^ 


(4.7) 


ai AlK 2 +a 2 A 2 7 K 2 = 



«2MlM3 = g 



(4.8) a 1 jjL l KA i -\-a 2 fji 2 KA 7 = 



a 2 jLt 3 KA 4 =- 


(4.9) 


a 1 x 2 r 4 +« 2 x s r 7 =i+| 


(4.10) 


« 1 # 2 r 4 +a 2 2Pr 7 =i+^ 


(4.11) 



Equations (4.5) are the usual equations for the integration parameters of a third order 
Runge-Kutta method. From (4.9) it follows that KA±=ii 2 and hence by (4.6) KA 7 =ix 2 . 
From (4.10) and (4.11) it is seen that KA l =j. Equations (4.6) through (4.11) may then 
be replaced by 



M2 

~K 



Mi 



a 1 r4+a 2 r 7 =^2 (l+3j) 



(4.12) 



No attempt will be made to discuss all possible solutions of (4.5) and (4.12). However, 
a few obvious facts will be pointed out. From (4.5) it is seen that mi^0> M3^0> a 2 7 £ 0. Like- 
wise Ti^O. It is clear that A lt A 4 , A 7 , r 4 and r 7 are independent. Since the matrix 

n o o o"\ 

10 
J) ai a 2> 
is of rank three any set of third order Runge-Kutta integration parameters may be used. 
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One could at this point make an analysis of various types of simple solutions of (4.12). 
However, in the next section it will be seen that there are some other considerations. Two 
sets of parameters are given for comparison purposes. 



(4.13) 



<*0, «1, «2, Hi, 


ix 2 , us (arbitrary solutions of (4.5)) 






, Mi 


*7 — TS — Ag 






►■ 


> l+3j 
As 6a 2yi K* 


"2 = ^3 = X5= 


=x 6 =x 9 =o 




> 






«o=0 


3 
«i=4 




1 

«s=4 




2 
M.=3 


M2=0 




M3=l 




x-- 7 ' 


x _2 
4 3K 


-x, 


x _ 2 






X2=X 3 =X 6 = 


=x 7 = 


= Xg=X 9 : 


= 0. 



(i+3i) 



(4.14) 



5. The Truncation Error 

In this section the truncated terms of the various expansions will be considered in some 
detail. A term shall be said to be T(j a K^h y ) if it is of the iorm.j a K^h?j where / is a function in- 
dependent of J, K and h. It would be desirable for all of the fourth order terms truncated in 
the integration of (1.2) to be T(A 4 ). The procedure would be pointless if some of the truncated 
terms are T(K*h*). It will be seen that there are terms which are T(Kh*) and T(j 2 K~ l W) which 
cannot be simultaneously eliminated by any choice of extrapolation parameters. 

All of the truncated terms from y mK +j(t mK+j +i) —y mKi . j+ i which are not T(A 4 ) are listed below. 
It is assumed that ki=jlK, A i =fx 1 /K, A 7 =jjl 2 /K. The terms from y mK+j (t mK+j+ i), d { (j) and 
d 2 (j) are grouped in that order. 

B 4^ (^+g^+g,«)^+ (1+8i t 12g,r ' ) OM 



24 



24 



(l + 4j) r r F , (1 + 4J+6J 2 ) (1 + 4J-6J 2 ) . 



» l K*T 4 (G xx F+G xy G+G xt )F 1 +IPhG s F x F 1 +^ 

f x 2 K 2 T 7 (G xx F+G xy G+G xt )F 1 +K^ 7 G x F x F 1 + ^K 3 A 7 G x F 2 + 

The original statement was that y(t) did not depend strongly on x(t) or varied more rapidly 
than x(t). This statement may be replaced by the following explicit assumptions on F and G: 
Assumption 1: The following inequalities are valid: 

\F,\<\Q,\IK, p=x,y,l,2. 

This assumption implies that y(t) does vary more rapidly than x(f). 
Assumption 2: G x is T(Kr l ). 
It is seen from (4.12) that T 4 and r 7 are T(jK~ 2 ). 



156 



Willi assumption 1 and the preceding remark it is seen that there are two groups of terms 
which are not T(h A ). They are 



(l+4j+6f) 1 

24 hx 2 ' 2 



G X F 2 , ~K^G X F 2 , ^IPA&Fz 



and 



(1+4J-6J 2 ) 
24 



G X F V G U K\K<t> i -jT i )G x F y G 1 , K?{K<t> t -3Y 1 )Q x Ffi i . 



These terms would be eliminated if the following equations held : 

1+4J+6J 2 A , . 
12g3 =a 1 A i +a 2 A 7 

1± ^£=a i (K<t>t-jr i )+a 2 (K*-,-T 7 ). 

Tlicse two equations may be transformed into 

l+4j+6,f _ 



UK' 



a 1 A4-|-«2^7- 



l + 8j+6,f 

24^3 — =«i<P4+a2<P7- 



(5.1) 
(5.2) 



Unfortunately, it is not possible to satisfy (4.12), (5.1) and (5.2) simultaneously. The 
equations 

.l+3j 



air 4 +a 2 r 7 = 



6K 2 

1+4J+6J 2 
12K 3 

1 + 8J+6J 2 



24iP 



are incompatible. It is possible to satisfy (4.12) and either (5.1) or (5.2). Since the terms 
involved by (5.1) and (5.2) are T^ErW) and T(Kh*), respectively, the smallest truncation error 
results when (4.12) and (5.2) are satisfied. 

If (5.2) or (5.1) are to be satisfied, then some of the simplicity of the extrapolation coeffi- 
cients is lost. Two sets of parameters are given below; the first satisfies (5.1) and the 
second satisfies (5.2). 



Mi, M2> Ti, 72 ? 73 (arbitrary Runge-Kutta integration parameters) 



_1 



Mi 



Xi — ^> X 4 — fjr — X5— \6 



_ l + 3j X 6T2 



K " 4 K " 5 " 6 V 6a iyi K 2 1 



_ l+4j+6i 2 -2 7l g(l + 3i) _ M 2 x , _ x _ x _ n 

X6_ 12^,72(72-7,) '~K X 2 -X 3 -X 8 -X 9 -0 

Mu M2 f 7i, 72, 7s (arbitrary Runge-Kutta integration parameters) ' 



x _.7 \ — a_x x 

Al — -jt A4 — TV- — A 6 — A 6 



= l+3j X 67 2 
iaaiK 2 7j 



x , l + 8j+6j» & . , _. _. _ fl 

X6 ~24Z 3 a 1 7 1 73 7 ~K X 2 -X 3 -X 8 -X 9 -0 



(5.3) 



(5.4) 
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These equations are somewhat more complicated than those given in section 4. 

With assumption 2 it is seen that there are three groups of terms which are not T(h*). 
They are 



24 



G X F X F U K^G X F X F U K^ n G x F x F Y 



and the two groups found with assumption 1. These terms would be eliminated if the following 
equations held. 

1+4J+6J 2 



2AJS? 

l+8j+6j 2 _l + 8j+12g 2 r 1 



O^+C^Ay 

^104+^207. 



24K 3 24iP 

Hence if I\ is set equal j 2 /2K 2 the same equations are found as with assumption 1. 

6. The Second Approach 

In this section the extrapolation parameters will be determined so as to make the trunca- 
tion error of the extrapolation as small as possible. This procedure is similar to the analysis 
for the Runge-Kutta integration of one equation. 

The extrapolation of x(t) to any value t mK to any value t mK -\-r is given by 

XmK+ Xo 0) fco + Xj (r) ki +X 2 (t) k 2 

where k , k t and k 2 are from (2.1). X (r), Xi(r) and X 2 (r) are determined by equating the co- 
efficients of h in the expansion of the error equal to zero. The resulting equations are 

Xo(r)+X 1 (r)+X 2 (r) = r/^/i (6.1) 

\{r) yi +\ 2 {r)r2=\r 2 l{KhY (6.2) 

Xi(r)7?+X 2 (T)7 2 2 =|r 3 /(^) 3 (6.3) 

Mr)ym=l^/(Khy. (6.4) 

Equation (6.1) results from equating the coefficients of h to zero, (6.2) results from equating 
the coefficients of h 2 to zero, and (6.3) and (6.4) result from equating the coefficients of h 3 to 
zero. In general only three of these equations can be satisfied One would naturally choose a 
solution which satisfies both (6.1) and (6.2). 

These equations may be written for the extrapolation parameters. The resulting equations 
are: 

Al= K Ai=> E A7= f (6 - 5) 

r _ i 2 r _ o"+mi) 2 r r __ (j+^y r ,„ fi x 

~2K* 4— " 2K 2 2K 2 ~ 

A - ? A . O'+Mi) 8 A A - Q'+^) 3 A (R 71 

Al -3K* Ai ~^K* Al A7 ~^K^~ Al (67) 
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Equations (6.6) and 6.5) imply that (4.12) is satisfied. The converse is not true. It is easily 
shown that (6.7) implies 



and that (6.8) implies 



A , A 6f+4j+4fa M i+<* 2 /4) 

a 1 A i +a 2 A 7 = 12g3 



,, . 6i 2 +4i+4(a lM ?+^) 



Therefore (6.7) implies that the G X F 2 term is T{j 2 K~ l h^) and (6.8) implies that the G x F y G l tennis 
TiKh*). Hence (6.7) and (6.8) have the same effect on the truncation error as (5.1) and 
(5.2) although (6.7) and (6.8) do not actually imply (5.1) and (5.2), respectively. 

It is rather surprising that T 1 =j 2 /2K 2 does not appear among the equations derived in sec- 
tions 4 and 5 with assumption 1. It is certainly plausible to include this equation in any set 
of equations taken to determine the extrapolation parameters. 

Two sets of parameters are given, the first satisfies (6.5), (6.6), and (6.7) and the second 
satisfies (6.5), (6.6), and (6.8). 



Mi, M2, 7i, 72, 73 (arbitrary Runge-Kutta integration parameters) 



Xi — 7?"~^ 2_ ^ 3 > 



x 4 — -jr— x 5 — x 6 



X7 jr ^8 \) 



\ 2 



_ 2f-Zy 2 Kf 
"6if 3 7i (7i—7 2 ) 



. _ 2Q'+Mi) 3 -37 2 g(i+Mi) 2 , 
6^ 3 y 1 ( 7l -72) 2 ' 

_ 2(,7+ M2 ) 3 -37 2 g(i+M2) 2 , 

6K 3 7i(7i-7 2 ) 



Xu= 



= 2f-ZKf 7l 

~6K 3 y 2 (y 2 —y l ) 

2(7+Mi) 3 -3g(j+Mi) 2 T 1 
6K 3 72 (72-7i) 

2(j+ M2 ) 3 -3^(j+M2) 2 7i 



6iP72(T 2 -7i) 



y (6.9) 



Mb M2, 7i, 72, 73 (arbitrary Runge-Kutta integration parameters) 



Xi — -jr— X 2 — X 3 , 



X4 — js — X 5 — X 6 , 



*7 — TV- — A 8 — A 9 , 



x 2 = 



Jhf 



~2K 2 y l (iK 3 y\y 3 



, _ Jj+M,) 2 7 2 (j+M,) 3 , 
As 2K 2 7 t " 6^ 3 7?Y3 2 ' 

x _ (j+M 2 ) 2 7 2 (j+M2) 3 . 
As_ 2K% " 6K 3 yh 3 2 ' 



3 0K 3 7,7 3 

x __(i+M,) 3 x 
> _ (i+M 2 ) 3 . 

X9 -6^77r 3 ~ 3 



(6.10) 



If only (6.5) and (6.6) are to be satisfied there is considerable simplification in the extrapola- 
tion coefficients. The following set satisfies (6.5) and (6.6): 



Xi=-£-X 2 



M2 

~K 



2K*y, 



> (i+M2) 2 , 
X8 -^^t7 _X2 



x 3 =o 



X - Ml -X X - O'+^l) 2 > 5,_n 

A 4 — 7^ — A5 A5— lyjSZ A 2 A6— O 



X 9 =0 



(6.11) 
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7. Fourth^Order Formulas 



An analysis similar to that of sections 4, 5, and 6 may be made for a fourth order integration 
method. Only the results of such an analysis are given here. The equations for a fourth 
order integration method are : 



#0 — Kht [X m K, VmK) W) 

k l =KhF(x mK +y 1 k 0) y mK +yiho, t m K+yiKh) 

k 2 =KhF(x mK +y 3 k 1 + (y 2 —y d )k , y mK +ydh 1 + (y 2 —y 3 )h , tmK+l2Kh) 

h=KhF(x mK +y 6 k 2 +y 5 k 1 + (y 4 —y 5 —y 6 )k , y m K+y&h 2 +y 5 h 1 + (y i —y 5 —y 5 )h 0) t mK +yiKh) 

h l =KhG(x mK +y 1 k , VmK+yih, tmK+yiKh) 

h 2 =KhG(x m K+ysk 1 +(y 2 —y 3 )k 0j y m K+yzh l -\-{y 2 —y 3 )h {)l t mK +y 2 Kh) 
h 3 =KhG(x mK +y & k 2 +y 5 k 1 + (74— y 5 — y Q )k , y m K+yQh 2 +y 5 h 1 + (74—75—76)^0, t nK +y A Kh) 

X(m+1)K — X m K ~l A)#0 I &# 1 "1 P2#2 + P3"?3 



XmK+j — >£_/ \i(j)ki-i-rX„, 



do(j) — h(*\X m K+j,ymK+j,tmK+j) 



di(j)=hG(z mK + J + y £} \iki-5,y m K+j+Hido,t mK +j+iiih) 



i=5 



d2(j)=hG(x mK +j+ y £] \iki- g ,y mK+j + n 3 d 1 + (ju 2 — ^)d ,t mK+j + ix 2 h) 



i=9 



dz(j)=hG(x mK+j + ^2 ^A-13,2/ot2S:+^ + M6^2 + M5^1+(M4— M5~ ^)d 0) t mK+j + fxjl) 



y mK+j+1 = y mK +j + a o d + <*i^i + «2^2 + oc z d z 



The integration parameters must satisfy: 



«0+«l + «2 + a3=l 
«lMl + «2M2 + «3M4 = I 
«1 Ml + «2M2 + «3M4 = 3 
«2MlM3+ «3M2M6+ ^3^^= i 



«lMl + «2M2 + «3/4 = I 

«2MlM2M3+ «3M2M4M6+ «3MlM4M5= : 8 

«2MlM3+«3M2M6+«3MlM5 = lV 



«3MlM3M6=2T 
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(7.1) 



The extrapolation parameters must satisfy the following equations in order for the 
integration method to be fourth order: 



a 1= =^ of 1 r 5 +a 2 r 9 +a 3 ri3=-j7^/ 



where 



~K 



a 1 A 5 +a 2 A9+a 3 A 1 3= 



1+4,7+6/ 
12K* 



M2 3+8? 

= K a i/ x ir 5 +a2M2r 9 +a3M4ri3=24^2 

\ ^ *.l ^ . ^ 1 + 8J+6J 2 

l3 ~*"/? «105 + «209 + «39l3 = 24/T3 

j 2 1 + 4J 

r i=2^2 «2M3r 5 +a3M5r 5 +«3M6r9==24^2" 



y (7.2) 



A n =X w +X w +i + X„+2 + X n + 3 A n = 7fX n+ i+7l\ n+ 2 + 7!Xn+3 

r n = 7l^»+l +72X n+ 2+74X w+ 3 0n = 7l73X n + 2 +(7l75 + 7275)^n+3 

The truncated terms have not been examined as in section 5 due to the extremely tedious nature 
of such an examination. 

The equations analogous to (6.5), (6.6), (6.7), and (6.8) are: 



Aj- 



Ti 



K 



-J— 
2K 2 



~K 



A 9 = 



M2 

K 



(i+Mi) 2 

2K 2 



-Ti r 



A -JL A - O'+Pi) 3 

n i — ory-3 n s— QZ ^ 3 -^1 ^9 



.(j+th) 



M4 

"if 



_ 0'+M4) 2 r 

1 13— o TS'l V 1 



2K* 



3K 3 



0,=7 



05 



3K 3 
.(j+fhl 



—<t>l 09 



3iJT 8 

.(J+M2) 3 



a, a 13 =(4±^ 3 - 



"6K 3 * 6_ 6K 3 ^ ™~ 6ii: 3 

The relation between the X^ and the A„, T n , A re and 4> n is 



01 0i: 



62P 



(7.3) 
(7.4) 
(7.5) 
(7.6) 



A, 




1 


1 


1 


1 




x» 


r„ 







7i 


72 


74 




^n+1 


A, 







7? 


7i 


7? 




\+2 


<t>n 










7i73 


(7i75+7 2 7e) 




^n + 3 



This relationship can be inverted and 



X,i 




A w 


^n + 1 
X w +2 


=^1 


An 


\*+3 




<t>n 
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where A is the matrix with elements a^lp^ 

P = 7i[(Y2 — Ti)72(Ti Y5 + Y2 7c) +Ti73 74(Ti — 74)] 
Gn = P «i2 = (7? — 72) (7i 75 + 72 7e) — 7i73 (ti— 7*) Gi8= (72 — 7i) (7i75 + 72 7e) + 7i7s (71 — 74) 

«i4= (71 — 72) (72—74) (71 — 74) 



a 2 i = ^22 = 72(7175 + 72 76)— 7l73 7l ^23 = 7173 74 — 72(7175 + 72 76) 



<h\ = a 32 = — 7i (7i 7s + 72 Ye) 



^33 = 7i(7i75 + 72 7e) 



a 41 = a 42 = 7i73 ^43=— 7i73 

If Runge-Kutta-Gill integration parameters are used, then 



^24 = 72 7 4 ~ 7 2 74 
^34 = 7l74(7l— 74) 
«44 = 7l 72(72 — 7l). 



A-- 



1 





10 



0.58578646 
3.414214 
-1 



2.8284272 
-6.8284272 
2 





-6.8284272 
6.8284272 




8. Experiments 



Experiments have been made with these formulas on three pairs of differential equations. 
They are 



dx . 
dt= x/2 

^=xcos(25t) 



dx 
di 



=y° 



dy_ 13x 

dt~a b (t+\) 12 

dx 



3(0) = 1 

2/(0) = 1/1250.5 

z(0)=a 6 /13 
2/(0)=a 



(8.1) 



(8.2) 



dt 

dy_ 

dt~ 



1 + (1.5X10" 4 ) cos (1500 z(0)=0 



=xcos (250 



y(0) = .00016. 



(8.3) 



Equation (8.1) is a typical example of a system suited for split Runge-Kutta. Equation 
(8.2) is an example of a system not suited for split Runge-Kutta. Yet even in this situation 
there may be some advantage to using it if the solution for y(t) is needed with much more 
accuracy than that for x(t). In (8.3), y(f) is much more rapidly varying than x(t), yet the 
small oscillation in x(t) gives it a comparatively large second derivative. The experiments for 
each of these equations are discussed in considerably more detail in sections 9, 10, and 11. 

No experiments were made with parameters derived from assumption 2 in section 5. 
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9. First Example The Typical Case 



The solution of (8.1) is 



x{t)=e 



— fi t/2 



(9.1) 



y(t) = (}{ cos 25J+25 sin 25t)e t/2 /625.25. 



(9.2) 



There are three basic variables in the solution of (8.1). They are (i) the method of solution, 
(ii) R, (iii) h. The methods of solution are numbered 1 through 6 corresponding to equations 
(4.13), (5.3), (5.4), (6.11), (6.9), and (6.10), respectively. The same integration parameters 
were used throughout. They were: 



a =A)=2/9, ai = ft=l/3, a 2 =ft=4/9, 7] = Ml = l/2, 72 = M2 =3/4, 73= 



: M3 = 



= 3/4. 



In figure 1 the errors in the y integration are studied as a function of K for various methods 
of solution. The integration interval h is fixed at .01 and the error is the maximum error for 
0</<l. In general this occurs close to f=l. The value of the error of the x integration 
at /=1 is also given. Here the variation of K is equivalent to varying the integration step. 
This error behaves like K 2 - 7 in the range of values considered. 



1 1 r 




Figure 1. Maximum integration errors as 
a function of K and the extrapolation 
method. 



v-. 
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There is very little difference in the results for the first three methods and only one curve 
is given. These methods proved to be inferior by a large factor compared with the second 
three methods. 

There is little difference in the results for the fourth and fifth methods. Only the results 
for the fourth method are plotted. This lack of difference can be explained. Method 4 
corresponds to satisfying eqs (6.1) and (6.2), and method 5 corresponds to satisfying (6.1), 
(6.2), and (6.3). Now (6.3) arises from setting the coefficient of F 2 equal to zero and for this 
system F=ax and F 2 =0. Hence there is no improvement when method 5 is used rather 
than method 4. 

On the other hand method 6, which satisfies (6.1), (6.2), and (6.4) is a definite improvement 
over method 4. Equation (6.4) arises from setting the coefficient of F x Fi-\-F y G\ = e t/2 /8 equal 
to zero and hence an improvement would be expected. 

There is another phenomenon apparent in figure 1 and a very significant one. That is the 
fact that the error does not decrease with K after a certain point. Indeed, the same accuracy 
for y(t) may be obtained with Kh=0.35 as with Kh==0.01 if method 6 is used. This phenom- 
enon is explained by the fact that there are two factors that contribute to the error in the 
y(f) integration. One of them is the error in extrapolation of x(t) and the other is the error 
inherent in the actual integration of y(t). For some values of K the extrapolation errors are 
dominant and for other values the integration errors are dominant. When the extrapolation 
errors are no longer significant then improvement can be made only by decreasing h. 

In figure 2 the errors are studied as a function of h with K fixed. This was done only for 




Figure 2. Maximum integration errors as 
a function of h and the extrapolation 
method. 
K is fixed at 25. 



.008 



.012 
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the first and fourth methods. Figure 2 is self-explanatory except to point out that things are 
not always as one thinks they should be in numerical analysis. There is no particular explana- 
tion for the fact that the error does not decrease monotonically with h, it just happens that way. 
Figures 3 through 7 are plots of error curves for various situations. Figure 3 is for method 1 , 
K=25 and h=.01. Figure 4 is for method 4, K=25 and h = .01. Figure 5 is for method 4, 
K=10, h = .01. It is identical with the curves with K=5, 1. Figure 6 is for method 6, R=25, 
h=.0l and it is identical with the curves with K=10, 5, 1 and very similar to the curve with 
K=35. Figure 7 is for method 1, K=--25, h=.01/\ / 2 ~ .0084. 



20XIO -8 




Figure 3. Error curve for method 1 with 
K=25 and h = 0.01. 



10X10 



5X10 




Figure 4. Error curve for method 4 with 
K=25 and h = 0.01. 




Figure 5. Error curve for method 4 with 

K=10andh = 0.01. 

Identical curves are obtained for smaller values of K. 




Figure 6. Error curve for method 6 with 
K=25 and h= 0.01. 
Identical curves are obtaired for smaller values of K 
and for if =35 the error curve is very similar, 




Figure 7. Error curve for method 1 with 

~K=25 andhtxO.0084. 

It is seen that the relatively larger error here is due to 

the chance shape of the curve and not to any basic lower 

accuracy. 
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10. Second Example — An Unprofitable Case 

The solution of (8.2) is 

x (t)=--a 6 (t+iyyi3 

y(t) = a(t+l)\ 



(10.1) 
(10.2; 



All of the results discussed below are for a=.l. 

In figure 8 the errors are studied as a function of K for various methods of solution. 



io- 4 - 



T 



1 

■ 1 ERRORS IN 
///X(t)ATt=.2 



'/,' 



V 



V 



/ 

/ / 

METHOD 4 + 5 > "V ERRORS IN Y (t) 

/ / / AT t =.2 

METHOD 6^/ V 

/' ' 

INCORRECT // / 
METHOD 6 ^' 



METHOD 4 + 5 
METHOD 6 




Figure 8. Integration error at t = 0.2 as a 
function of K and extrapolation method. 
It is seen that the y(t) error depends directly on the 
x(t) error. The maximum error accurs at £=0.2. 



25 



Methods 1, 2, and 3 are not used. Methods 4 and 5 give very similar results while method 6 
does somethat better than either of them. Also included is an "incorrect method 6" which 
gives the best results of all for larger K's. In this method X 8 was computed from 



„-!* 



-M2) 2 , i 72(j + Mi) 3 
2K 2 7l 6KVl7z 

instead of its true value given in (6.10). 

The fact that method 6 is better than methods 4 and 5 cannot be explained as previously. 
Indeed, by that argument method 5 should be the best and method 4 the worst with method 6 
in between. Kecall that method 5 satisfies (6.3) which is obtained from setting the coefficient 
of F 2 equal to zero, likewise method 6 is obtained from setting the coefficient of F x Fi+F y 6i 
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equal to zero. Now in this case F 2 =120a 6 (t+l) 10 and F x F l + F y G l =6a (i (t+1) 10 and hence 
method 5 should be the best. This points out that the truncation error is not always a good 
measure of the integration error. 

It is seen that the integration error does not level off as K decreases. This is due to the 
fact that the major portion of the error in the y(t) integration is due to the errors in the inte- 
gration of x(t), as differentiated from the errors in the extrapolation of x(t). This is borne out 
by figures 9 and 13. In figure 13 the y(t) integration error is plotted for the system of equations 



dx 
dt 



=a\t+l) 12 



^=rsx/a 5 (t+iy 2 

which have the same solution as (8.2). Method 6, K=20 and A=.002 is used. This leads to 
much more accurate values of x(t) which in turn lead to a much more accurate y(t) integration 
even though the extrapolation accuracy has not improved. 

In figure 9 the error curves for methods 4, 5, and 6 are compared with K=20, h= 0.002 
and for 0<t<.2. 

The error curve for method 4 with R=l, A = 0.002 is given in figure 10. 

In figure 11 two error curves are given which have the same interval of integration for 
the x(t) equation. Both use method 6, one with K=l and h=.02 and one with K=10 and 




Figure 9. Error curves for methods 4, 5 
and 6 with K = 20 and h = 0.002. 

The improvement of method seems to lie in the fact 
that the error remained positive much longer than for 
the other two methods. 



24XI0" 9 



Figure 10. Error curve for method 4 with 
K=l and \\ = G.002. 

This corresponds to the normal Runge-Kutta integra- 
tion procedure. 




12X10 



h=.02 k=l 




Figure 11. Error curves for method 6 with 
K and h taken to be (10, 0.002) and 
(1, 0.02). 

These curves emphasize the fact that the accuracy is 
determined by Kh and not h alone. 



22X10 




Figure 12. The unusual error curve for the 
incorrect method 6 with K= 10 and h = 
0.002. 

The improved answers are apparently due to the fact 
that the error remains positive for a long period. 



552370—60- 
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Figure 13. Error curve of the modified 
equation for method 6 with ~K=20 and 
\i = 0.002. 
The modification removed the dependence of the x{t) 

equation on y(t). The maximum error is decreased by a 

factor of 10. 



^=.002. The close similarity of the two curves again emphasize that the accuracy of the 
x(t) integration determines the accuracy of the y(t) integration. 

The unusual error curve for the "incorrect method 6" is given in figure 12 with K=10 
and A = 0.002. 
p> Some studies were also made with a=2 but the overall behavior was the same. 



1 1 . Third Example — Noise 



The solution to (8.3) is 



x(t)=t /10 + 10" 6 sin 150 
fy(0=4X10" 8 [sin 25/-6 sin 3 25/+(48/5) sin 5 25f— (32/7) sin 7 25*] +0.004 [/sin 25/+.04 cos 25/]. 

This system is very similar to (8.1) with the addition of a rapid oscillation on x(t). However 
the experimental results bear no resemblance to those for (8.1). In fact they seem to be com- 
pletely chaotic. Although the errors in the y(t) integration are all relatively small and fairly 
uniform in size, there appears to be no correlation between the maximum error and if or method. 
An examination of the error curves reveals the same lack of trend. Even the errors in the 
x(t) integration do not behave rationally as a function of K, which is equivalent to the inte- 
gration interval. Indeed the largest error occurs for i?=4. 

In figure 14 the y(t) integration errors are plotted as functions of K and method. The 
x(f) integration errors, which are independent of the extrapolation method, are also plotted 
as a function of K. 

Figures 15 to 18 are error curves for various situations. 

12. Numerical and Efficiency Considerations 

The ingredients of an ideal situation for these formulas are as follows: (1) F(x,y,t) is com- 
plicated and difficult to evaluate, (2) the solution for y(t) is relatively insensitive to errors in 
x(t); i.e., the main error source is the inherent inaccuracy of the y(t) integration. The first- 
ingredient magnifies the efficiency and the second is necessary to attain high accuracy. The 
second ingredient may be replaced by "the accuracy desired for y(t) is lower than that required 
for x(t)." It should be remembered that numerical integration is a delicate business in general 
and these formulas are particularly so and hence each case should be examined on its own 
merits. 

A comparison will now be made between the work required for split Kunge-Kutta and 
the normal Runge-Kutta methods. 

Unless the parameter K is changed very often the extrapolation parameters should be 
computed only once. For the integration of (1.1) and (1.2) from t to to+Rh the following 
tables give the required number of operations for normal Runge-Kutta and split Runge-Kutta 
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Figure 14. Maximum integration errors as 
a function of K and extrapolation method. 

The lines only join thedafa points and are not intended 
to represent the actual behavior of the curves between 

data points. 




Figure 16. Error curve for method 4 with 
K = 6 and h = 0.01. 

The error oscillates with approximately constant 
amplitude and a period similar to that of figures 5, 15. 
and 17. 



I8XI0" 8 




•160 xiO" 




Figure 17. Error curve for method T 4 with 
K = 5 and h = 0.01. 
The error oscillates as in figure 16 except it is now 
oscillating about the line; error = (1.2^+0.2) 10~«, instead 
of zero. 



Figure 15. Error curve for method 6 with 
K = 2 and h = 0.01. 
This curve was reproduced for methods 5 and 6 with 
K=l and ft =0.01. Apparently the majpr portion of the 
error oscillates with a period of about 0.23. Upon this 
then there is superimposed an error with much smaller 
magnitude and a much shorter period. Compare the 
period of the error with figure 5. 



QO x 10' 




Figure 18. Error curve for method 4 with 
K=4 and h = 0.01. 
The same oscillation period is present and the oscilla- 
tion is about zero, but the error curve appears to be 
"amplitude modulated." 
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methods. Table 1 is for third order integration and table 2 is for fourth order integration. 
From these tables it is seen that if a large K may be used then there is a large saving in com- 
putation even if F is extremely simple. If F is complicated then the efficiency is even greater. 



Table 1 





Add. 


Mult. 


deval- 
uation 


G eval- 
uation 


Normal Runge- 
Kutta.. 


22K 


2SK 


3K 


SK 






Split Runge-Kutta 
method 4 


14CK-+D 


17if-HQ 


3 


SK+2 


Split Runge-Kutta 
methods 5, 6 


1GK+U 


19CET+1) 


3 


3K+2 



Table 2 





Add. 


Mult. 


JP eval- 
uation 


G eval- 
uation 


Normal Runge- 
Kutta.. . __ 


SSK 


4QK 


\K 


\K 






Split Runge- 
Kutta ___ _ 


2hK+21 


291^+34 


4 


4iT+3 







(Paper 64B3-32) 
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